load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/substring.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/WRF_dyy_contributed.ncl"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;            User Config                   ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

if (.not. isvar("configFile")) then
    configFile = "/public/users/wrf/Codes/WRF_scripts/WRF_Post/StationLists/sswdn_stnlist.csv"
    print("configFile is undefined, using default: " + configFile)
else
    print("Using configFile: " + configFile)
end if

lines = readAsciiTable(configFile,1,"string",1)
delim = ","

stnnum    = ndtooned(str_get_field(lines,1,delim))
stnlist    = ndtooned(str_get_field(lines,2,delim))
stnlon    = ndtooned(tofloat(str_get_field(lines,3,delim)))
stnlat    = ndtooned(tofloat(str_get_field(lines,4,delim)))
stndomain = ndtooned(str_get_field(lines,5,delim))
nstns     = dimsizes(stnlist)

; default output time interval is 30min
if (.not. isvar("output_interval")) then
    output_interval = 0.25 ; 15 min
    print("output_interval is undefined, using default = " + output_interval)
else
    print("Using output_interval = " + output_interval)
end if

; test if input_file is given
if (.not. isvar("input_file")) then
    input_file = "/public/users/wrf/wrfout/2012103120/wrfout_d02_2012-10-31_12:00:00"
    print("input_file is undefined, using default = " + input_file)
else
    print("Using input_file: : " + input_file)
end if

filename = systemfunc("basename " + input_file)

; test if output_prefix is given
if (.not. isvar("output_prefix ")) then
    output_prefix = "~/wrfpost/test/"
    print("output_prefix is undefined, using default = " + output_prefix)
else
    print("Using output_prefix: " + output_prefix)
end if

if (.not. isfilepresent(output_prefix)) then
    system("mkdir " + output_prefix)
end if

; check start_date
if (.not. isvar("start_date")) then
    start_date = substring(filename, 11, 23)
    print("start_date is undefined, using default = " + start_date)
else
    print("Using start_date: " + start_date)
end if

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;           Open NetCDF File               ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
a = addfile(input_file + ".nc", "r")

domain = substring(filename, 7, 9)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;              Process Time                ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
times = a->Times
hours = WRF_Times2Udunits_c(times, 0)
wrf_dt = hours(1) - hours(0)
; UTC -> CST
hours = hours + 8 ; 
timeStr_cst = ut_string(hours, "%Y-%N-%D_%H:%M:%S")
ntimes = dimsizes(hours)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;               Read SSWDN                 ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
sswdn = wrf_user_getvar(a,"SSWDN",-1)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;          Loop for each station           ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

do istn = 0, nstns - 1, 1
    ; if not in the desired domain
    if (.not. isStrSubset(filename, stndomain(istn))) then
        continue
    end if

    ; Setting Output Filename
    output_filename = output_prefix + "/test_cubic_spline_" + stnnum(istn) + "_" + stnlist(istn) +  "_" + domain +  "_" + ut_string(hours(0), "%Y-%N-%D_%H") + ".txt"

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;             Lon/Lat to XY                ;
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    res = True
    res@returnInt = False   ; return real values
    loc = wrf_user_ll_to_ij(a, stnlon(istn), stnlat(istn), res) - 1

    sxc0 = floattointeger(loc(0))
    syc0 = floattointeger(loc(1))

    syc1 = syc0 + 1
    sxc1 = sxc0 + 1

    yy = (/syc0, syc1/)
    xx = (/sxc0, sxc1/)

    print ("X/Y location is: " + loc)

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ; crop station area from whole data array          ;
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    sswdn_crop = sswdn(:, syc0:syc1, sxc0:sxc1)
    
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ; calculate interpolated stational data            ;
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    sswdn_station = ndtooned(linint2(xx, yy, sswdn_crop, False, loc(0), loc(1), 0))

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;       do the cubic spline interpolation          ;
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    knots       = ntimes
    new_timenum = fspan(hours(0), hours(ntimes-1), tointeger((hours(ntimes-1) - hours(0))/output_interval + 1))
    copy_VarAtts(hours, new_timenum)
    printVarSummary(hours)
    new_sswdn   = dyy_interp_ts(hours, sswdn_station, new_timenum)
;    new_sswdn   = csa1(hours, sswdn_station, knots, new_timenum)
    printArray1  = timeStr_cst + ", " + tostring(sswdn_station)
    printArray2  = ut_string(new_timenum, "%Y-%N-%D_%H:%M:%S") +  ", " + tostring(new_sswdn)
    printArray = array_append_record(printArray1, printArray2, 0)
    print(printArray)
    asciiwrite(output_filename, printArray)

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ; output precipitation tendency for dt hours       ;
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;    l_cnt = 0
    ; write header line
;    line_array(l_cnt) = "Station Number, Station Name, Time, Shortwave Radiation"
;    l_cnt = l_cnt + 1
;    do it = dt, ntimes-1, dt   ; TIME LOOP
;        line = stnnum(istn) +", "+ stnlist(istn) +", " + timeStr_cst(it) + ", " + sprintf("%5.2f", sswdn_station(it))
;        line_array(l_cnt) = line
;        l_cnt = l_cnt + 1
;    end do
    ;asciiwrite(output_filename, line_array)
    ;print("Creating: " + output_filename)
end do ; loop for station
